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Abstract —Since many decades, there is a general perception 
in literature that the Fourier methods are not suitable for 
the analysis of nonlinear and nonstationary data. In this pa¬ 
per, we propose a Fourier Decomposition Method (FDM) and 
demonstrate its efficacy for the analysis of nonlinear (i.e. data 
generated by nonlinear systems) and nonstationary time series. 
The proposed FDM decomposes any data into a small number 
of ‘Fourier intrinsic band functions’ (FIBFs). The FDM presents 
a generalized Fourier expansion with variable amplitudes and 
frequencies of a time series by the Fourier method itself. We 
propose an idea of zero-phase filter bank based multivariate FDM 
(MFDM) algorithm, for the analysis of multivariate nonlinear 
and nonstationary time series, from the FDM. We also present an 
algorithm to obtain cutoff frequencies for MFDM. The MFDM al¬ 
gorithm is generating finite number of band limited multivariate 
FIBFs (MFIBFs). The MFDM preserves some intrinsic physical 
properties of the multivariate data, such as scale alignment, trend 
and instantaneous frequency. The proposed methods produce the 
results in a time-frequency-energy distribution that reveal the 
intrinsic structures of a data. Simulations have been carried out 
and comparison is made with the Empirical Mode Decomposition 
(EMD) methods in the analysis of various simulated as well 
as real life time series, and results show that the proposed 
methods are powerful tools for analyzing and obtaining the time- 
frequency-energy representation of any data. 

Index Terms —The Fourier decomposition method (FDM); 
Fourier intrinsic band functions (FIBFs) and analytic FIBFs 
(AFIBFs); zero-phase filter bank based multivariate FDM 
(MFDM); Empirical Mode Decomposition (EMD). 

1. Introduction 

HE time-frequency representation (TER) is a well es¬ 
tablished powerful tool for the analysis of time series 
signals. There exist many types of time-frequency (TE) anal¬ 
ysis methods, e.g. linear (the short-time Eourier transform), 
quadratic (the Wigner-Vihe distribution) and Wavelet trans¬ 
forms. The TER is achieved by formulation often referred 
as time-frequency distribution (TED) and provides insight 
into the complex structure of a signal consisting of several 
components. These approaches have many useful applications, 
however the analysis of nonstationary signals are not well 
presented by these methods. 

Recently developed Empirical Mode Decomposition 
(EMD) d has provided a general method for examining 
the TED, and has been applied to ah kinds of data. The 
EMD is an adaptive signal analysis algorithm for the 
analysis of nonstationary and nonlinear signals (i.e. signals 
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generated from nonlinear systems). The EMD has become 
an established method for the signal analysis in various 
applications, e.g. medical studies El-ll, meteorology HI, 
geophysical studies ini and image analysis 0. The EMD 
decomposes any given data into a finite number of narrow 
band intrinsic mode functions (IMEs) which are derived 
directly from the data, whereas other signal decomposition 
techniques (like the Eourier, Wavelets, etc.) incorporate 
predefined fixed basis for signal modeling and analysis. The 
Ensemble EMD (EEMD) is a noise-assisted data analysis 
method developed in IH to overcome the timescale separation 
problem of EMD. The Multivariate EMD (MEMD) developed 
in a is a generalization of the EMD for multichannel data 
analysis. The Compact EMD (CEMD) algorithm is proposed 
in na to reduce mode mixing, end effect, and detrend 
uncertainty present in EMD and to reduce computation 
complexity of EEMD as well. The IMEs, generated by EMD, 
are dependent on distribution of local extrema of signal 
and the type of spline used for upper and lower envelope 
interpolation. The traditional EMD uses cubic spline for 
upper and lower envelope interpolation. The EMD algorithm, 
proposed in ifTTll to reduce mode mixing and detrend 
uncertainty, uses nonpolynomial cubic spline interpolation to 
obtain upper and lower envelopes, and have shown in ca 
that it improves orthogonality among IMEs. 

The energy preserving property is important for any kind 
of transformation, and it is obtained by the orthogonal de¬ 
composition of signal in various transforms like the Eourier, 
Wavelet, Eourier-Bessel, etc. The energy preserving property 
is especially important for the accurate and faithful analysis 
of three dimensional time-frequency-energy distribution of a 
signal. The EMD algorithms, proposed in |[T3ll . ensure orthog¬ 
onality or energy preserving property or both in decomposition 
of signal into IMEs and refereed to as energy preserving EMD 
(EPEMD). 

In spite of considerable success of EMD, ah of the EMD 
algorithms are based on empirical, heuristic and ad hoc proce¬ 
dure that make them hard to analyze mathematically, and EMD 
may suffer from mode mixing, detrend uncertainty, aliasing 
and end effect artefacts There is a lack of mathematical 
understanding of the EMD algorithms, e.g. IMEs dependence 
on the number of sifting and the stopping criteria, convergence 
property and stability to noise perturbation. In spite of ah 
these limitations, EMD is the widely used nonstationary data 
analysis method. Therefore, in this paper, EMD is used as a 
reference to establish the validity, reliability and calibration of 
the proposed method. 

Since many decades, there is an understanding in the litera- 
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ture (e.g. ID, ca, lED) that the Fourier methods are, directly, 
not suitable for nonlinear and nonstationary data analysis, and 
various reasons (e.g. linearity, periodicity or stationarity) are 
provided to support it. The Fourier transform is valid under 
very general Dirichlet conditions (i.e. signal is absolutely 
integrable with finite number of maxima and minima, and 
finite number of finite discontinuities in any finite interval) thus 
include nonlinear and nonstationary signals as well. Therefore, 
in this study, we explore and provide algorithms to analyze 
nonlinear and nonstationary data by the Fourier method termed 
as the Fourier Decomposition Method (FDM), which generates 
small number of Fourier intrinsic band functions (FIBFs). It 
is already well established that the Fourier method is best 
tool for spectrum analysis and, in this study, we show that 
the Fourier method is also a best tool for time-frequency 
analysis and processing of any signal. The power of the Fourier 
transform can also be realized from the fact that the analytic 
representation and, hence, the Hilbert transform of a signal 
are, inherently, present in the Fourier transform. 

In this paper, we also propose a method, which captures the 
features of the MEMD, using a zero-phase filter bank (ZPFB) 
approach to construct the multivariate FIBFs and residue com¬ 
ponents. This multivariate FDM (MFDM) algorithm generates 
matched multivariate FIBFs and residue through zero-phase 
filtering. Thus, we propose an adaptive, data-driven, ZPFB 
based time-frequency analysis method. 

For the adaptive data analysis approach, the most difficult 
challenge has been to establish a general adaptive decomposi¬ 
tion method of analysis without a priori basis. In this study, we 
propose the FDM and MFDM general adaptive data analysis 
methods that are inspired by the EMD methods and their filter 
bank properties lITSl . ifT^ . This paper is organized as follows: 
In section II the EMD algorithm is briefiy presented. We 
propose the Fourier decomposition method (FDM) in section 
III. We propose the ZPFB based multivariate FDM algorithm 
in section IV. Simulation results are presented in section V. 
Finally conclusions are presented in section VI. 

II. Brief overview of the EMD algorithm 

There are various tools for nonstationary data processing 
such as the spectrogram; the wavelet analysis; the Wigner- 
Ville distribution; evolutionary spectrum ifTvll : the empirical 
orthogonal function expansion (EOF) (or principal component 
analysis or singular value decomposition); Synchrosqueezed 
wavelet transforms dSl; the EMD, etc. 

The EMD is a well established signal decomposition method 
that decomposes nonlinear and nonstationary data into a set 
of finite band-limited IMEs and residue through the sifting 
process. The decomposed signal x{t) is expressed as the sum 
of i IME components plus the final residue as 

£ £-\-l 

x{t) = '^yi{t)+re{t) = '^yi{t), (1) 

i=l i=l 

where yi{t) is the IME and ri{t) = is fi¬ 

nal residue. The IMEs admit amplitude-frequency modulated 
(AM-EM) representation |[T^ (i.e. yi{t) ~ ai(t) cos(0i(t)), 
with ai{t), > 0 Vt) and well-behaved Hilbert 


transforms m. Eor any IME yi{t), its Hilbert transform yi{t) 
is defined as convolution of yi{t) and l/irt, i.e. yi(t) = 
^ dr and the Hilbert transform emphasizes the local 

properties of yi{t). An analytic signal Zi{t) can be repre¬ 
sented by Zi{t) = yi{t) ^ jyi{t) = ai(t) exp (j0i(f)), where 
O'iit) = and <pi{t) = tan-^[yi(t)/yi{t)] are 

instantaneous amplitude and phase of yi{t). The instantaneous 
frequency (IE) of yi{t) is defined as: uJi{t) = 0^(f) = 
Vi The physical meaning of IE uJi{t) con¬ 

strains that (t)i{t) must be a mono-component function of 
time. The Bedrosian and Nuttall theorems CSl, EOl fur- 
ther impose non-overlapping spectra constraints on the pair 
of a signal yi{t) = ai(f) cos(</>i(f)). 

All IMEs must satisfy two basic conditions: (1) In the 
complete range of time series, the number of extrema (i.e. 
maxima and minima) and the number of zero crossings are 
equal or differ at most by one. (2) At any point of time in 
the complete range of time series, the average of the values 
of upper and lower envelopes, obtained by the interpolation of 
local maxima and the local minima, is zero. The first condition 
ensure that IMEs are narrow band signals and the second 
condition is necessary to ensure that the IE does not have 
redundant fluctuations because of asymmetric waveforms |T| . 

HI. The Eourier Decomposition Method 

We propose a class of functions, termed as the Eourier 
intrinsic band functions (EIBEs), belonging to C°^[a,6], here 
with the following formal definition. 

Definition 1: The Eourier intrinsic band functions (EIBEs), 
yi{t) G C°^[a,6], are functions that satisfy the following 
conditions: 

(1) The EIBEs are zero mean functions, i.e. yfit) dt = 0. 

(2) The EIBEs are orthogonal functions, i.e. yi{t)yj{t) dt = 
0, for i 7 ^ j. 

(3) The EIBEs provide analytic EIBEs (AEIBEs) with instanta¬ 

neous frequency (IE) and amplitude always greater than zero, 
i.e. yi{t)+jyi{t) = withai(t), > 0, 

yt. 

Thus, the AEIBEs are monocomponent signals and, physically, 
the IE has meaning only for monocomponent signals, i.e., 
signal has only one frequency or a narrow range of frequencies 
varying as a function of time 1^ . Thus, the EIBE is sum of 
zero mean sinusoidal functions of consecutive frequency band. 

The main objective of this study is to obtain unique repre¬ 
sentation of multicomponent signal as a sum of constant and 
monocomponent signals, i.e. signals which can be represented 
by the following model 1241 : 

M 

xit) = Z] 2/* w+(2) 

i=l 

where n{t) is a noise representing any residue (constant or 
trend) components and the yfit) are M single component 
nonstationary signals, which in our proposed framework would 
be the EIBEs, defined above. 

The necessary conditions d, for a basis to represent a 
nonlinear and nonstationary time series, are completeness. 
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orthogonality, locality, and adaptiveness. The FIBFs, intrin¬ 
sically, follow all the necessary conditions by virtue of the 
decomposition. 

The available data are usually of finite duration, nonsta¬ 
tionary and generated from the systems that are generally 
nonlinear. Let x{t) be a time limited + Tq] real valued 
signal which follows the Dirichlet conditions. We construct 
the periodic signal as XTo(f) = — sTq) such that 

x{t) = XTo(t)u;(t), where w{t) = 1, for ti < t < ti + Tq 
and zero otherwise. The Fourier series expansion of Xj^ (t) is 
given by 

oo 

XTo(t) = ao + '^[dk cos{kujot) + bk sm{kuJot)], (3) 

k=l 

where frequency (rad/s) ccq = ^ XTo(t)df, 

a/c = XjQ{t) cos{kujot) dt and bk = 

^ Xroit) sm{kuJot) dt. We write (O as 

^ oo 

Xioi-t) = «0 + 2 (jktoot) + 4exp {-jkuot)], (4) 

^ fe=i 

where Ck = {dk — jbk) and = {ak jbk). From dH), it is 
clear that 

XTo(t) = ao + Re{zj^{t)}, (5) 

where analytic function 

oo 

%.(*) = E Ck exp (jkuJot) (6) 

k=l 

is complex conjugate of Zioit) - E^iCfcexp {-jkuiot) and 
Re{zjQ{t)} denotes the real part of ZjQ{t). We write ZjQ{t) as 

M 

^To W = Ll (7) 

i=l 

where, in forward search (i.e. low to high frequency scan) 
of AFIBFs, ai(f)exp(ji'(/)i(f)) = X]fe=i exp (jfcwof), 
a2(f)exp(j</>2(f)) = Er=(iVi+i) Cfcexp(ifcwof), 

aM{t) exp (jcpMit)) = E^(jvm-i+i) i.e. 

Ni 

ai{t)exp{j(t)i{t)) = Cfeexp(jA:wof), (8) 

k = Ni—i-\-l 

with A^o = 0 and Nm = oo. The FIBFs are the real part 
of AFIBFs presented in Eq. ([U). In order to obtain minimum 
number of AFIBFs in low to high frequency scan (LTH-FS), 
for each i, start with (A^i_i + 1) and append more term till 
we reach the maximum value of Ni such that (A^i-i + 1) < 
Ni < oo and 

ai{t),LOi{t) = —> 0,Vf- (9) 

It is easy to observe that such a decomposition is always 
possible. 

Similarly, in reverse search (i.e. form high to 
low frequency scan (HTL-FS)) of AFIBFs, we 
obtain ai{t)exp{j(f>i{t)) = Cfc exp (jfcwof), 

a2(f)exp(j</>2(f)) = Efe=iV2^^Cfeexp(jfcwof), ■■■, 


aM(f)exp(j(?!>M(f)) = E1=T ' ^Efeexp(ji'fca;of), and 

the lower and upper limits of sum in Eq. (Eli would change to 
k = Ni to {Ni-i — 1), respectively, with Nq = oo, Nm = 1. 
Here, we start with (A^i_i — 1), decrease and select minimum 
value of Ni such that I < Ni < {Ni-i — 1) and Eq. (|9]) is 
satisfied for i = 1, • • • , M. 

Observe that (|7]) has precisely the form that in the lit¬ 
erature m is termed as a generalized Fourier expansion. 
Moreover, it is worth noting that this representation is com¬ 
plete, orthogonal, local, adaptive and purely Fourier based. 
Thus, we have obtained a generalized Fourier expansion of 
a time series in Eq. (|7]) by the Fourier method itself. The 
variable amplitude and the IF have improved the efficiency 
of the expansion by expanding the signal into finite number 
of analytic FIBFs, in (|7]), and enabled the expansion to 
accommodate nonstationary data. Thus, we have obtained a 
variable amplitude and frequency representation, whereas, the 
classical Fourier expansion provides the constant amplitude 
and fixed-frequency representation. 

For each FIBFs, the amplitude ai(t) and IF fi{t) are 
functions of time, therefore, we define the three dimensional 
{t, /i(t), di{t)} time-frequency distribution of amplitude as the 
Fourier-Hilbert spectrum (FHS) H{f, t). The marginal Hilbert 
spectrum (MHS), derived from Hilbert spectrum, is defined 
in im. Similarly, here we define the marginal Fourier-Hilbert 
spectrum (MFHS) from the FHS as follows: 

h{f)= r H{fR)dt. (10) 

Jo 

The marginal Fourier-Hilbert spectrum offers a measure of 
total amplitude (or energy) contribution from each value of 
frequency in a probabilistic sense. The frequency in either 
H{fR) or h{f) has a different meaning from the Fourier 
spectral analysis fii . The presence of energy at each frequency 
in MFHS h{f) means that, in the total duration of the signal, 
there is a higher likelihood for such a wave (FIBF) to have 
appeared locally. The frequency in the MFHS indicates only 
the likelihood that an oscillation with such a frequency exists. 
The exact occurrence time of that oscillation is given in the full 
Fourier-Hilbert spectrum. We can also define the instantaneous 
energy density, which can be used to measure the fiuctuation 
of energy with time, as 

rfM 

E{t)= / ( 11 ) 

Jo 

where /m is a maximum frequency of signal. From Eq. (O, 
we obtain the energy of signal x{t) (or power of signal Xj^ (t)) 
by the ParsevaTs theorem sls Ex = | [^n + ^n] 

from Eq. ^ energy of the analytic signal (or power of signal 
Zto(^)) as Ez = + ^n]’ therefore, relation between 

Ex and E^ is given by 

Ex = dQ . (12) 

Hence, energy of zero mean signal is half of the energy of its 
analytic signal. 

Since in practice the continuous time signals are, generally, 
discretized for further processing by a computing device, so 
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we present the FDM for discrete signal. Let, x[n], be a discrete 
signal of length N. Using the discrete Fourier transform 
(DFT), we can write x[n] as 

= L] (13) 

/c=0 

where X[k] = X x[n] exp( ) is the DFT of 

signal x[n]. Let N be an even number (we can proceed in 
the similar fashion when N is an odd number), then X[0] and 
X[y] are real numbers; and we can write x[n] as 


N-l 

x[n] = X[0] + ^ X[k] exp{-^^^) + X[y] exp(j7rn) 
k=l 

+ (14) 

/c=^ + l 

Since x[n\ is real, therefore, zi[n\ = X[k] 

is complex conjugate of Z 2 [n] = X[/c] 

and we can write m as 

x[n] = X[0] + 2Re{z^[n]} + X[y](-1)”, (15) 

where Re{zi[n]} denote the real part of zi[n]. Now, we write 
analytic signal zi[n] as 

i2iTkn A 

X[fc]exp(:C_) = Yai[n]exp{j4>i[n]), (16) 

k=l i=l 

where, in forward search (i.e. low to high frequency 

scan) of AFIBFs, we obtain ai [n] exp [n]) = 

Ef=i^[fc]exp(^), a2[n]exp{jMn]) 
E^=(jVi+i)^Wexp('^^^)’ aM[n]exp{j(f>M[n]) = 

EfcL(jvM-i+i) N )’ 

ai[n]exp{j4)i[n\= Y ^[fc] exp( '^^^^^ ), (17) 

k=Ni-i + l 


with A^o = 0 and Nm = {^ ~^)Xn order to obtain minimum 
number of FIBFs in LTH-FS, for each i, we scan from 
(Ni-i +1) to (Y — 1), obtain maximum value of Ni such that 
{Ni_i-\-l) < Ni < (X — ^) and phase 0i[n] is a monotonically 
increasing function, i.e. uji[n] = (0i[n 1] — (piln]) > 0 or 


„,[„] = ( ■».[■■+ ll-0.[«-l| ) J 


Vn 


(18) 


and ai[n] > 0 for i = 1, • • • , M and Vn. Observe that such a 
decomposition always exists. 

Similarly, in HTL-FS for FIBFs, we obtain 

ai [n] exp [n]) = E& m exp(^), 

a2[n]exp{jMn]) = ES ■■■, 

aM[n]exp{j(f>M[n]) = Efe=T~'~^^ ^ 

the lower and upper limits of sum in Eq. dnii will change to 
k = Ni to {Ni_i — 1), respectively, with Nq = ]\fj^ = 

Here, for each i, we scan from (A^i-i — 1) to 1, obtain the 


minimum value of Ni such that l<A^i<(A^^_i — 1) and 
phase (pi [n] is a monotonically increasing function. 

Thus, FDM provides two views, low to high frequency and 
high to low frequency view, of the signal and generate two 
set of time-frequency-energy distribution. Depending on the 
signal, both view may be same or sometimes they reveal 
two different types of features of the signal. The FDM is 
summarized in Algorithms A and B. The FDM with Fourier 
transform (FT) and discrete time Fourier transform (DTFT) is 
summarized in appendix. 

Algorithm A: The FDM algorithm (LTH-FS) to obtain 
AFIBFs, for 2 = 1, • • • , M with A'o = 0 and Nm = 

STEP 1. Obtain X[k] = FFT{x[n\}. 

STEP 2. Set AFIBFi = E^=(iv,_,+i) X[k] expi^^) = 
ai[n]exp{jpi[n]), obtain maximum value of Ni such 
that {Ni-i + 1) < Ni < {X — and phase pi[n] 
of AFIBFi is a monotonically increasing function, i.e. 

= Vn. _ 

Algorithm B: The FDM algorithm (HTL-FS) to obtain 
AFIBFs, for 2 = 1, • • • , M with Nq = X and Nm = 1- 
STEP 1. Obtain X[k] = FFT{x[n]}. 

STEP 2. Set AFIBFi = X[k]exp{i^^) = 

ai[n]exp{jpi[n]), obtain minimum value of Ni such that 
^ ^ Ni < {Ni_i — 1) and phase pi[n] of AFIBFi 
is a monotonically increasing function, i.e. uji[n] = 

^ ^ ^ Q 

IV. Multivariate Eourier Decomposition Method 

Prom (El and (113, we observe that the operation that 
generates the PIBPs is nothing but the Eourier based zero- 
phase filtering (ZPP). This is another motivation, in addition 
to the PB properties of IMPs, to use the Eourier or other 
methods of zero-phase filtering to decompose any data into 
a set of PIBPs. The ZPP of a real valued signal x[n] by 
zero-phase filter (hz[n]) can be obtained by two methods: 
(1) Convolution method, i.e. y[n] = x[n] * {hz[n]) ^ 
Y[k] = X[k]Hz[k], where hz[n] = h[n] ^ h[—n] and 
IIz[k] = |Ff[/c]p, where h[n] is a real sequence. (2) The 
Pourier method, i.e. set IIz[k] = 1 at desired frequency 
band and Hz [k] = 0 otherwise, obtain output by the inverse 
DPT, i.e. y[n] = X[k]Hz[k]ex.p{j27rkn/N), where 

X[k] = exp{-j2Trkn/N). 

We use ZPP that does not shift the essential features of the 
signal, and propose a multivariate PDM (MPDM) algorithm to 
generate multivariate PIBPs (MPIBPs) and residue as follows: 
Apply zero-phase high pass filtering (ZP-HPP) with cutoff 
frequency fd to each of the components of the P-variate (P- 
channel) time series {xp{t)}p^i and obtain first set of MPIBP 
ypi{t). The first set of residue is obtained as follows: 

rpi{t)=Xp{t)-ypi{t), p=l,2,---,P. (19) 

Apply ZP-HPP with cutoff frequency fc 2 to set of residue 
Vpi (t) and obtain second set of MPIBP yp 2 (t). The second set 
of residue is obtained as 


Tp2 (^) — Tpi (t) yp2 (t) p — 1, 2, • • • ^ P. 


( 20 ) 
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We can repeat this ZP-HPF procedure ^ times and obtain final 
set of MFIBF ypi{t) and residue (with cutoff frequency fd) 

rpi{t) = - ypi{t) p=l,2,---,P. (21) 

Through the addition of (1^ and (1^ we obtain expres- 
sion, similar to o, for P-variate time series as 

£ 

Xp{t) = '^ypi{t)+rpe{t), p=l,2,---,P. (22) 

i=l 

When we use the Fourier based zero-phase filtering, as in 
Eq. (Hi, to obtain MFIBFs, first two conditions of FIBFs are 
fully satisfied and the third one is approximately satisfied 
(i.e. satisfied in all practical sense), obviously, it can not 
be guaranteed simultaneously to all P-channel data. This is 
similar to MEMD algorithm problem where in derivation of 
multivariate IMFs, first condition of IMF is not imposed (91 . 

The question is, how to obtain cutoff frequencies (CFs) 
/cl 5 /c 2 , • • • 5 fci corresponding to zero-phase high pass filters 
hz 2 {t)^ • • • , hz£{t)l There is lot of fiexibility and are 
various ways to select CFs, e.g. dyadic (i.e. fd = fc2 = 
^• •) fc£ = ^ , where /m is the maximum frequency of 
a signal x{t) and for the sampled signal, maximum frequency 
is (^) half of the sampling frequency), non-dyadic, uniform 
and non-uniform CFs. We can take the Fourier transform of 
signal x{t) to obtain its spectrum details and make strategy to 
decide CFs. 

For narrowband signal, we define ratio of center frequency 
(fc) to bandwidth (BW) as 

m = fa/if Hi - fLi), fa = {fni + fLi)/^^, (23) 

where fni is the highest frequency and fa is the lowest 
frequency of ith band of a filter bank. From ([2^ we obtain 

fa = [(2m - l)/(2m + l)]fHi, rn > 1/2. (24) 

From ([23]) and (l24b . we observe that the ratios, for the 
consecutive ith and (i + l)th bands, of center frequencies fa, 
CFs fci and BWs {fni — fhi) can be taken as a constant, i.e. 

fci/fci+i = fci/fci+1 = , = 1. (25) 

[jHi+1 - jLi+l) 

From (23]), (|24]) and (25]), we obtain I = (2m + l)/(2m — 1) 
or m = (l/2)(/ + !)/(/ — 1) with I > 1, and as [m ^ oc, 
/ ^ 1], [m ^ 1/2, / ^ oc]. Here, we have liberty to select 
any suitable value of / or m, and greater the value of m (or 
lesser the value of 1) narrower the band, whereas in the case 
of dyadic FB / = 2 and m = 1.5 are fixed values. If required, 
we can vary the value of m (or /) for each band rather than 
taking the fixed value. Thus, we here propose the compact 
and elegant way to decide CFs as summarized in Algorithm C. 
Algorithm C: An algorithm to obtain cutoff frequencies fd. 

1. Select suitable value of m and set fni = Fs/2. 

2. Set fci = [(2m - l)/(2m + 1)]/^,,. 

3. Set /iTi+i = fci. 

4. Repeat step 2 to 3 for i = 1, 2, • • • , 

In MFDM, we can use zero-phase low pass filtering (ZP- 

LPF) in place of ZP-HPF to decompose signal in order of 
residue to first MFIBFs, i.e. rp^(^), ^p^(t), • • • ,^pi(t). We 


Xj(t) + nj(t) X2{t) +11 ^( 1 ) Xj(t) + n3(t) x^W + n/t) 



Fig. I: MFDM applied to a quadri-variate tone-noise mixture, generating perfectly aligned 
intrinsic modes in all the four channels. 


use zero-phase filtering as it preserves salient features (e.g. 
maxima, minima, etc.) in the filtered time waveform exactly at 
the time where those features occur in the unfiltered waveform, 
whereas conventional (non zero-phase) filtering shifts the 
features in the signal and hence cannot be used. The zero- 
phase filtering of time series can be obtained through the 
finite impulse response (FIR) or infinite impulse response (HR) 
filters. 

Similar to the MEMD and noise-assisted MEMD (NA- 
MEMO) ca, this MFDM algorithm produces the equal 
number of scale-aligned MFIBFs for all channels and pre¬ 
serving joint channel properties that make it suitable for 
direct multichannel modelling. The FDM does not suffer from 
mode mixing, detrend uncertainty and end effect artefacts as 
extraction of FIBFs does not depend on distribution of local 
extrema across the range of signal. 

Table jl] presents comparisons among Fourier, Wavelet, 
EMD 123 and FDM methods in data analysis. 

V. Simulation results 

The online available MATLAB software of MEMD (2T]|, 
EMD and EEMD fT2\ have been used in simulation results. 


A. Multivariate data decomposition 

We used quadri-variate time series signal, which is sum¬ 
mation of sinusoids (with combination of frequencies /i = 
4:Hz,f2 = SHz^fs = 16Hz,f4, = 32Hz) and Gaussian 
white noise of zero mean and standard deviation of 0.2., i.e. 
Xjit) = Eie[i, 4 ] + nj{t) (for j = I,-- - ,4). A 

32Hz sinusoid is present to first, second and fourth channels; 
a 16Hz sinusoid is present to first, second and third channels; 
a SHz sinusoid is common to all channels; a AHz sinusoid 
is present to first, third and fourth channels. On the same 
machine, computation time for MEDM is 0.45 sec. and for 
MEMD is 69.5 sec. in this simulation. The MEDM algorithm, 
similar to MEMD, generating perfectly aligned intrinsic bands, 
as shown in Eigure[T] in all the four channels, whereas, MEDM 
is computationally more efficient. 
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TABLE I: Comparisons among the Eourier, Wavelet, EMD-Hilbert and proposed EDM methods in data analysis. 



Fourier 

Wavelet 

EMD-Hilbert 

FDM 

Basis 

a priori 

a priori 

adaptive 

a priori ^ adaptive 

Frequency 

convolution: 

global 

convolution: 

regional 

differentiation: 

local 

differentiation: 

local 

Uncertainty 

yes 

yes 

no 

no 

Presentation 

frequency- 

energy 

frequency- 

time-energy 

frequency- 

time-energy 

frequency- 

time-energy 

Nonlinear 

no 

no 

yes 

yes 

Nonstationary 

no 

yes 

yes 

yes 

Harmonics 

yes 

yes 

no 

yes ^ no 

Theoretical 

base 

complete 

complete 

empirical 

complete 




Eig. 2: Decompositions of a sinusoid mixed with intermittent interference through 
algorithm of (b) EDM (c) EMD (d) EEMD. 


ting in EMD algorithm. These issues are mitigated by EEMD 
and NA-MEMD. The decompositions of a signal x{t) through 
EDM algorithm is shown in Eigure Ob) without end effect 
artefacts. The MEDM is able to localize the mono-component 
sinusoid within a single EIRE and outperforming EMD Oc) 
and EEMDO^). On the same machine, computation time for 
MEDM, EMD and EEMD are 0.56 sec., 0.21 sec. and 77.18 
sec., respectively. The ensemble size for EEMD was N = 500 
with the 16.94 dB signal-to-noise power ratio. 

C. Time-Frequency-Energy Analysis 

Eigure [3] shows time-frequency-energy (TEE) estimates for 
a nonstationary signal mixture of a linear chirp and frequency 
modulated (EM) sinusoid, obtained using the EDM and EMD. 
There is a enhanced TEE tracking when using EDM with low 
to high frequency scan and other one (HTL-ES) is similar to 
plot obtained by EMD algorithm. 

D. Intrawave frequency modulation 

Eirst, we decompose the following signal that has intrawave 
frequency modulation and it is considered challenging because 
the instantaneous frequency itself has very high frequency 
modulation ll^ 

1 cos(327rt H-0.2 cos(647rt)) 

1.2 + cos{27Tt) ^ 1.5 + sin(27rt) 

Second, we consider a model wave 

x{t) = cos{ujt + esinccf) (27) 

that satisfies the following highly nonlinear differential equa¬ 
tion d 

— ^ ^ +[c(;+ecc cos(co’t)]^x(t) —[ecc^ sin(ci;t)] y^l — x‘^{t) = 0 
^ (28) 
with uj = 1 and e = 0.5. We demonstrate that, Eigure O [6] 
and [71 our method well applies to these challenging cases with 
good accuracy. These examples clearly demonstrate that the 
EDM can indeed analyze nonlinear signals and it is a nonlinear 
decomposition method. 


B. Intermittency and mode mixing 

The intermittency, in the time series, is a main cause of 
mode mixing [e.g. signal x{f) in Eigure Oa)] and mode split- 


E. Analysis of a white Gaussian noise 

Eigure [8] shows the TEE analysis of a white Gaussian noise 
(with zero mean, unit variance, 1024 samples and sampling 
frequency Eg = 100 Hz) obtained from the EDM and EMD 
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FDM: Time-Frequncy-Energy estimate of FIBFs (LowToHighFrequencyScan) 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Time (s) 

FDM: Time-Frequncy-Energy estimate of FIBFs (HighToLowFrequencyScan) 
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EMD: Time-Frequncy-Energy estimate of IMFs 


Fig. 3: The TEE analysis of nonstationary signals which is mixture of linear chirp and 
EM sinusoid: FDM (top and middle) and EMD (bottom). 



■^0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Time (s) 


Fig. 4: The signal x{t) (solid line), sum of DC and lowest frequency FIBF (dashed line) 
of (26). 



Fig. 5: The DC, FIBFl, FIBF2 and highest frequency component by FDM. The DC 
and highest frequency component correspond to term X[0] and X[^]( —1)”^ of Eq. 
O, respectively. Sum of DC, FIBFl, FIBF2 and highest frequency component exactly 
synthesize signal x{t) given by Eq. {26) . 


FDM: Time-Frequncy-Energy estimate of FIBFs (LowToHighFrequencyScan) 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Time (s) 



Fig. 6: The TEE analysis of signal x{t) given by Eq. (2^: FDM (top) and EMD (bottom). 


algorithm. Clearly, both LTH-FS and HTL-FS views of TFE 
is similar and complete data is decomposed in FIBFs. Fig¬ 
ure [9] shows the power spectral density (PSD) plot of same 
white Gaussian noise with the FDM and EMD algorithm. 
The FDM has dived the complete data in narrowband and 
orthogonal FIBFs. Both LTH-FS and HTL-FS views of PSD 
looking similar but FIBFs have different frequency band, e.g. 
approximately 40 Hz is cutoff frequency for one of the band 
in PSD (LTH-FS), whereas, it is mid frequency of the one 
band in other PSD (HTL-FS) view. There are enhanced TFE 
and PSD tracking when using FDM. 
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FDM: Time-Frequncy-Energy estimate of FIBFs (LowToHighFrequencyScan) 



0 5 10 15 20 25 30 

Time (s) 

EMD: Time-Frequncy-Energy estimate of IMFs 



Time (s) 



Fig. 7: The TEE analysis of signal x{t) given by Eq. Wt\ : FDM (top) and EMD (bottom). 

F. TFE Analysis of unit sample sequence 

The unit sample sequence is defined as (5[n — no] = 1 at 
n = riQ and zero otherwise. By using the relation z[n] = 
^ X{uj) exy){jujn) duj, we obtain the analytic representa¬ 
tion of x[n] = S[n — no] X{uj) = exp(—jccno) as 

. , sin(7r(n - no)) + j[l - cos(7r(n - no))] 

Z\n\ = -7-r- 

7r(n - no) (29) 

z[n] = a[n] exp(j0[n]), 

where real part of z[n] is 6[n — no] = ^ = 

I I’ ~ f ~ hence uj[n] = ^ which 

corresponds to half of the Nyquist frequency, i.e. ^ Hz. 
Figure [TOl shows the plots of real, imaginary part and absolute 
value of z[n] with no = 199, sampling frequency Fs = 100 
Hz, length N = 400. Theoretically, this clearly indicate that 
most of the energy of signal 6[n — no] is concentrated at time 
t = 1.99 sec. (no = 199) and frequency / = 25 Hz. Figure fT^ 
shows the plots of FIBFs obtained by the FDM and Figure [TT] 
shows the TFE analysis of unit sample sequence ^[n—no] from 
the FDM, EEMD and continuous wavelet transform (CWT) 
methods. This clearly indicate that the TFE plot obtained by 
the FDM method is same as theoretical estimation, whereas 
there is energy spread over a range of frequencies and lack of 
accuracy in TFE plot by the EEMD and CWT methods. 

The uncertainty principle is a consequence of the Fourier 
transform (or any other type of integral transform) pair in the 
frequency definition. Therefore, its limitation could only be 
applied to such integral transforms, in which time would be 


EMD: Time-Frequncy-Energy estimate of IMFs 



Time (s) 


Fig. 8: Tbe TFE analysis of tbe Gaussian white noise (with zero mean and unit variance): 
FDM (top and middle) and EMD (bottom). 


‘integrated out’ or smeared over the integral time limit. Con¬ 
sequently, if we avoid an integral transform in the frequency 
computation, time-frequency analysis of signal would not be 
bounded by the uncertainty principle 1^ . The IF is defined 
through differentiation of phase rather than integration and, 
hence, overcome the restriction of the uncertainty principle. 
One can obtain arbitrary precision on time and frequency reso¬ 
lution subject only to the sampling rate in case of discrete time 
signal. The uncertainty principle in signal analysis state that 
the finer the time resolution one wants, the cruder the resulting 
frequency resolution would be. However, this example clearly 
demonstrate that the FHS by FDM is indeed not limited by 
uncertainty principle and signal can be highly concentrated in 
both time and frequency domain. 
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FDM: Power spectral density of FIBFs (LTH-FS) 




FDM: Power spectral density of FIBFs (HTL-FS) 

501-^^^^^^^ 



_40ol -^^^^^^^^^- 

0 5 10 15 20 25 30 35 40 45 50 

Frequency (Hz) 


EMD: Power spectral density of IMFs 



Fig. 9: The PSD analysis of the Gaussian white noise (with zero mean and unit variance): 
FDM (top and middle) and EMD (bottom). 


G. Application to a earthquake signal analysis 

Earthquake time series signal is nonlinear and nonstationary 
data. The Elcentro Earthquake data (sampled at Fs = bOHz) 
has been taken from ll^ and is shown in Eigure [131 The 
critical frequency range that matter in the structural design is 
less than lOHz, and from the Eourier based power spectral 
density (PSD), marginal spectrum by EDM and marginal 
spectrum by EMD in Eigure [13] show that almost all the 
energy in this data is within lOHz. The instantaneous energy 
fluctuations by the EDM and EMD methods, as shown in 
Eigure [141 are similar in nature. The TEE distribution by the 
EDM, EMD and CWT methods are shown in Eigure [TS] and all 
three methods indicate maximum energy concentration around 
1.7Hz and 2 second. There is enhanced TEE tracking by EDM 
methods as it provide better details of how the different waves 


Fig. 10: The analytic representation of S[n — no] (with, no = 199, sampling frequency 
Fs = 100 Hz, length N = 400): Real part of z[n] (top). Imaginary part of z[n] 
(middle) and absolute value of z[n] (bottom). 


arrive from the epical center to the recording station, e.g. the 
compression waves of small amplitude but higher frequency 
range of 10 to 20Hz, the shear and surface waves of strongest 
amplitude and lower frequency range of below bHz which 
does most of the damage, and other body shear waves which 
are present over the full duration of the data span. 

H. Application to a speech signal analysis 

To illustrate the advantages of EDM in speech signal anal¬ 
ysis, we decompose a quasiperiodic voice signal (small-cap 
I vowel). In Eigure [Tb] the EIBEs 1/2 to seem to extract 
most of the noise, while to yj seem to express most of the 
signal information at high frequencies. Moreover, the EIBE yg 
captures accurately FO, the fundamental frequency (inverse 
of glottal pulse or pitch period length) of the voice signal. In 
contrast, EMD (Eigure [TT] (a)) presents mode mixing between 
modes two and three, between modes three and four, and 
between modes four and five. Also, the EMD is not able to 
catch Fq in any mode. In Eigure [T7](b), the EEMD’s is able to 
catch Fq in mode ^ 5 , although it can be observed that positive 
minima and negative maxima in the mode y^ and yj, there is 
clear violation to IME condition (1). Besides this problem, the 
major drawbacks in this case are huge computation complexity 
for this ensemble size and the reconstruction error that may 
be significant. The computation complexity can be reduced 
by selecting the smaller ensemble size but that increases 
the reconstruction error. Eigure [18] shows the plots of TEE 
obtained by the EDM, EMD and EEMD algorithms for same 
data (vowel ‘small-cap T). Clearly, there is enhanced TEE 
tracking when using EDM. The results in this application as 
well are in clear favor of the Eourier method proposed in this 
paper. 


VI. Conclusion 

In this paper, we have proposed: (1) The Eourier De¬ 
composition Method (EDM), for nonlinear and nonstation¬ 
ary time series analysis, which decomposes any data into a 
small number of ‘Eourier intrinsic band functions’ (EIBEs). 
The EDM is a generalized Eourier expansion with variable 
amplitudes and frequencies of a time series by the Eourier 
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FDM: Time-Frequncy-Energy estimate of FIBFs (FTH-FS) 
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EEMD: Time-Frequncy-Energy estimate of IMFs 
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Fig. II: The TEE analysis of unit sample sequence S[n — no] (with, no = 199, 
sampling frequency Fs = 100 Hz, length N = 400): FDM (top), EEMD (middle) and 
CWT (bottom). 


method itself. (2) The zero-phase filter bank based multivariate 
FDM (MFDM) algorithm, for the analysis of multivariate 
nonlinear and nonstationary time series, which is generating 
finite number of band limited multivariate FIBFs (MFIBFs). 
(3) An algorithm to obtain cutoff frequencies required in 
MFDM algorithm for zero-phase high or low pass filtering 
of multivariate signals. 

The fundamental and conceptual contributions of the this 
study are the Fourier based decomposition method (i.e. FDM) 
and the introduction of the FIBFs. The FIBFs form the basis 
of the decomposition that are complete, orthogonal, local 
and adaptive. The instantaneous frequencies of the FIBFs 
produce a time-frequency-energy distribution of any signal. 
A time-frequency-energy distribution of a signal is used in 
various fields of science and engineering for analysis of 


^ 2r 

3 

t 0 


c _2l-^^^^^^^- 

0 0.5 1 1.5 2 2.5 3 3.5 4 



Time (s) 

Fig. 12: The DC, FIBFl, FIBF2 and highest frequency component plots by the FDM of 
unit sample sequence 6[n — no] (with, no = 199, sampling frequency Fg = 100 Hz, 
length N = 400). 

physical phenomena and engineering systems. The proposed 
methods produce the final presentation of the results in a 
time-frequency-energy distribution that reveals the imbedded 
structures of a signal. Unlike the various EMD algorithms, the 
FDM and MFDM are mathematically well defined, supported 
by the well established theories of filter and the Fourier 
transforms. The FDM and MFDM methods do not suffer from 
mode mixing, detrend uncertainty and end effect artefacts 
as extraction of FIBFs does not depend on distribution of 
local extrema across the range of data. Simulation results 
demonstrate the efficacy of the proposed methods. 

Appendix 

A. The FDM for continuous time real function 

Let x{t) be a non-periodic, real function of time and follow 
the Dirichlet conditions, then the Fourier transform (FT) of 
x{t) is defined as 

/ oo 

x{t)eyi-p{-j2'Kft)(lt (30) 

-CX) 

and inverse Fourier transform is defined as 

/ oo 

X{f)ex.p{j2TTft)df (31) 

-OO 

It is easy to show that A(/) exp(jf27r/t) d/ = 

/~X(-/)exp(-j27r/i)d/. From Eq. ® -X(-f) = 

X*(f) and, hence, we rewrite Eq. dSTIl as 

nOO 

x{t)= / (/) exp(j27r/t) +X*{f) exp(-j27r/t)] d/ 

Jo 

(32) 

In this Eq., second term is complex conjugate of first term. 
As x{t) is a real function, we can write 

x{t) = Re{z{t)} (33) 

where analytic function z{t) = 2 X{f) ex.p{j2Fft) df 

and Re{z{t)} denote real part of function z{t). We write 
Eq.([3^ as 

pOO ^ 

2/ X{f)ex.p{j2Trft)df = 'Fai{t)ex.p{j4>i{t)) (34) 

i=l 
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Elcentro Earthquake May 18,1940. North-South Component E(t) by EDM 






Marginal spectrum by EDM 



Marginal spectrum by EMD 



Fig. 13: Plot of (top to bottom): (a) Elcentro Earthquake May 18, 1940 North-South 
Component, (b) Fourier based power spectral density (PSD), (c) Marginal spectrum by 
EDM, (d) Marginal spectrum by EMD. 

where (with /o = 0, /m = oo) 

rfi 

ai{t)exp{j(j)i{t)) = 2 X{f)exp{j2Trft)df, (35) 
Jfi-i 


Fig. 14: Plot of the instantaneous energy density E(t) (top to bottom) using the: (a) 
EDM and (b) EMD. 


for i = 1, • • • , M. To obtain minimum number of AFIBFs in 
low to high frequency scan, for each i, start with fi-i, increase 
and select the maximum value of fi such that fi-i<fi<oo 
and 

aiit),fiit) = ^ ( 36 ) 

Similarly, in high to low frequency scan, the lower and 
upper limits of integration in ([35]) will change to fi to fi-i, 
respectively, with /o = (X),/m = 0, and we can obtain 
minimum number of AFIBFs by selecting the minimum value 
of fi such that 0 < fi < fi-\ and Eq. (|36|) is satisfied. 

From Eq. (|32|) and ([33]), it is easy to obtain relationship 
between the energy of original signal x{t) and energy of its 
analytic signal z{t) as 

Ez 

Ex = 

that is, the energy of analytic signal is twice of the energy of 
original signal. 

The IF characterizes a local frequency behavior as a function 
of time. In a dual way, the local time behavior as a function of 
frequency is characterized by the group delay (GD) : ti{f) = 
— ^ . The GD measures the average time arrival of the 

frequency /. In general, the IF and GD define two different 
curves in the time-frequency plane. Similar to the IF process, 
for causal signal x{t), we obtain 

ai{f)exp{-j4>i{f))= f x{t)exp{-j2TTft)dt, (38) 

Jti-i 
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Fig. 15: Plot of TEE (top to bottom) using the: (a) EDM (b) EMD and (c) CWT. 


with to = 0,tM = oo such that ai{f),ti{f) = 
0, V/, for i = 1, • • • , M. 


B. The FDM for discrete time real function 


Fig. 16: Vowel ‘small-cap P sound with sampling frequency Fg = 16000 Hz (top left) 
and its Fourier spectrum (top right). The FDM generates highest frequency component 
(yi), FIBFs {y 2 to yn) and DC term (r/ 12 ). The fundamental frequency Fq is captured 
accurately in FIBF yg. 
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Fig. 17: IMFs (yi to ^ 12 ) generated by (a) EMD and (b) EEMD (with zero mean and 
0.2 standard deviation of the added noise of normal distribution and the ensemble size 
of 300) algorithms for vowel ‘small-cap T with Fg = 16000 Hz. The EMD is not able 
to catch Fq in any mode, and the EEMD is able to capture Fq accurately in IMF r/ 5 . 


Let x[n] be a non-periodic, real function of time and 
follow the Dirichlet conditions, then the discrete time Fourier 
transform (DTFT) of x[n] is defined as 

00 

^h = E x[n] exp(—jfccn) (39) 

— oo 

and inverse discrete time Fourier transform (IDTFT) is defined 
as 

1 f"" 

—'F~ / X{uj) ex.-p{jujn) duj (40) 

271 J — ^ 


It is easy to show that X(uj) exp(jujn) dec = 
/q X(—cc) exp(—jccn) dec. From Eq. X(—ec) = X*{ec) 
and, hence, we rewrite Eq. (001) as 


x\n\ = 


1 r r 

— [/ X{ec) exp{jecn) deed-f X*{ec) exp{—jecn) dec]. 
271 Jo Jo 


(41) 

In this Eq., second term is complex conjugate of first term. 
As x[n] is real function, we can write 


x[n] = Re{z[n]} 


(42) 
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FDM: Time-Frequncy-Energy estimate of FIBFs (HTL-FS) 



Time (s) 


function, i.e. 

ai[n],uji[n] = ((/)i[n + 1] - > 0, Vn. (45) 

Similarly, in high to low frequency scan, the lower and upper 
limits of integration in Eq. (l44l) will change to uji to cci_i, 
respectively, with ccq = 7r,co’M = 0, and we can obtain 
minimum number of AFIBFs by selecting the minimum value 
of uJi such that 0 <uJi < cc^-i and Eq. (Ell is satisfied. 
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Fig. 18: The TEE plot by FDM (top) EMD (middle) and (b) EEMD of same data (vowel 
‘small-cap T). There is enhanced TEE tracking when using FDM. 


where analytic signal z[n] = ^ X{uj) exj[){jujn) duj and 
Re{z[n]} denote real part of function z[n]. We write Eq. (l42l) 
as 


^ «7r M 

— / X(c(;) exp(jfci;n) dec = > ai[n] exp(jf0i[n]) (43) 


where (with ccq = 0, uom = tt) 

1 

ai[n]ex.p{j(l)i[n]) = - / X(cc) exp(jf'ccn) dcu, (44) 

^ JuJi-l 


for i = 1, • • • , M. In order to obtain minimum number of 
AFIBFs in low to high frequency scan, for each i, start with 
uji-i, increase and select the maximum value of uJi such that 
< tt and phase 0i[n] is a monotonically increasing 
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